function [y, yh, hx, mr, X] = sac(varargin)
%    [yan, sachist, sacx, mr] = sac(allsp, width, binw, delay, dur);
%
% Shuffled autocorrelation function
% Based on Louage et al., 2004
% X is an array of N responses x length(t) points in time. The times are in
% msec. Since the number of spikes in X may vary from trial to trial,
% X is assumed to be a cell array of spike times.
% The autocorrelation is calculated for every event in X that is not within
% is not within a single spike train.
% twin is the time window for the autocorrelation to be calculated (in
% msec)
% binw is the bin width for the histogram, in msec.
% delay is the delay to the start of measurement in X; in msec.
% dur is the duration of the measurement window in X; in msec.
% (these last two parameters allow you to pass a response and only analyze
% the portion that is relevant).
%
% y is the output autocorrelation function. It consists of a list of
% the correlation intervals, not corrected at all.
%
% yh is the histogram bin amplitudes, with bin width hx.
% This histogram is normalized.
% mr is the mean rate
%
% 10/5/04 Paul B. Manis, Ph.D.
% 10/5/04 some fixes for an spike trains.
% Note: calling the routine with letters as the first argument will generate the plots
% from Figure 2 of Louage et al.

global ALLCH DFILE


y = []; yh = []; hx = []; mr = 0; X={};

switch nargin
    case 0
        return;
    case 5
        X = varargin{1};
        twin = varargin{2};
        binw = varargin{3};
        delay = varargin{4};
        dur = varargin{5};



    case 1 % call as a test - we make up our data and plug it in, and at the end
        % we make our own plot. the first arg is the letter corresponding to
        % the figures in Louage et al., 2004, Figure 2.
        %
        cmd = varargin{1}; % capture command

        baseper = 4/3;
        stimdur = 1000;
        ntestrep = 20;
        switch(cmd)
            case 'A'
                for i = 1:ntestrep
                    X{i} = [baseper:baseper:stimdur];
                end;
                twin = 5;
                binw = 0.05;
                delay = 0;
                dur = stimdur;
                ddur = 10;
            case 'B'
                for i = 1:ntestrep
                    X{i} = [baseper:baseper:stimdur];
                    twin = 5;
                    X{i} = X{i} + 0.08*randn(size(X{i}));
                end;
                binw = 0.05;
                delay = 0;
                dur = stimdur;
                ddur = 10;
            case 'C'
                fprintf(1, 'There is no computation for C, only A,B,D,E,F and G\n');
                return;

            case 'D'
                for i = 1:ntestrep
                    X{i} = [baseper:baseper:stimdur];
                    X{i} = X{i} + 0.170*randn(size(X{i}));
                end;
                twin = 5;
                binw = 0.05;
                delay = 0;
                dur = stimdur;
                ddur = 10;
            case 'E'
                for i = 1:ntestrep
                    bl = [baseper:baseper:stimdur];
                    bp = randperm(length(bl));
                    X{i} = sort(bl(bp(1:floor(length(bp)/2)))); % elimnate half by random selection
                    X{i} = X{i} + 0.170*randn(size(X{i}));
                end;
                twin = 5;
                binw = 0.05;
                delay = 0;
                dur = stimdur;
                ddur = 10;
            case 'F'
                twin = 5;
                binw = 0.15;
                delay = 0;
                dur = stimdur;
                for i = 1:ntestrep
                    X{i} = stimdur*(sort(rand(120,1))); % corresponds to about 120 s/s
                end;
                ddur = 100;
            case 'G'
                twin = 5;
                binw = 0.15;
                delay = 0;
                dur = stimdur;
                sig = stimdur*(sort(rand(120,1)));
                for i = 1:ntestrep
                    X{i} = sig; % corresponds to about 120 s/s
                end;
                ddur = 100;

            otherwise
                return;
        end;


    otherwise
        return
end;


% now the SAC calculations.

lx = length(X);
% make sure input bin width and window are integer mulitples.
if(((twin/binw) - floor(twin/binw)) ~= 0)
    x=floor(twin/binw);
    twin = binw * x; % recalculate the window
end;
maxlx = 100;
maxn = 100000;
%fprintf(1, 'l = %d\n', lx);
if(lx > maxlx) % for testing purposes.
    lx = maxlx;
end;
yc=cell(lx,1); % Prellocation helps ALOT with speed
tic;
n = 1;
spcount = zeros(lx, 1);
for i = 1:lx
    xi = X{i}; % get spike train i
    kxi = find(xi >= delay & xi < delay+dur); % window the data to be analyzed
    xi = xi(kxi); % reduce data set
    spcount(i) = length(xi);
%    fprintf(1, 'i=%d  ns: %d n', i, spcount(i)); % let user know we're working
    tic;
    yj = NaN*zeros(maxn,1); % preallocate space for cross correlation
    n = 0;
    for j = 1:lx
        if(j ~= i)
            xj = X{j}; % get spike train j
            kxj = find(xj >= delay & xj < delay+dur);
            xj = xj(kxj); % reduce data set
            for ti = 1:length(xi) % cross correlate against all spikes in xi
                iti = xi(ti); % time of ti'th spike in spike train i
                xjk = find(xj >= iti-binw/2 & xj < iti+twin+binw/2); % find spikes in j, in the window relative to current spike
                % binw/2 are egde effect corrections for the histogram stuff below.
                nx = length(xjk);
                if(~isempty(xjk) & (n+nx) < maxn)
                    n = n+1;
                    yj(n:n+nx-1) =  [xj(xjk)-iti]; % compute the difference in time between spike in j and relative to current spike in i
                    n = n + nx-1;
                end;
            end;
        end;
    end;
    yc{i} = yj; % somehow this is faster to store this way...
%    fprintf(1, ' elapsed: %8.3fs  last n = %d\n', toc, n);
end;

ya=[yc{:}]; % turn cell array into nx1 linear array
y=ya(find(~isnan(ya))); % clean it up - and shorten.
if(length(y) == 0) % no spikes in this window
    return;
end;

% now calculate the normalized histogram.
% normalization is N*(N-1)*r^2*deltat*D
% where N is the number of presentations (lx), r is the rate (sp/sec),
% deltat is the bin width for the histogram (sec), and D is the duration of
% the trace (?).
% normalization goes to 1 as the time gets larger...
%
rate = sum(spcount)/(lx*dur/1000); % get average rate
%fprintf(1, 'Mean firing rate: %f\n', rate);
mr = rate;
nfac = lx*(lx-1)*rate*rate*(binw/1000)*(dur/1000); % correction factor
[yh] = histc(y, [0:binw:twin]); % generate histogram
hx = [0:binw:twin];
yh = yh/nfac; % to convert to rate, spikes/second
%m = mean(yh);
%fprintf(1, 'Mean value of yh: %8.3f, scale factor: %f\n', m, nfac);

if(nargin == 1) % if we are testing, generate a plot...
    hf = findobj('tag', 'sactest');
    if(isempty(hf) | (hf <= 0))
        hf = figure;
        set(hf, 'tag', 'sactest');
    end;
    figure(hf);
    clf;

    subplot(2,1,1);
    hold on;
    for j = 1:length(X)
        k = find(X{j} < ddur);
        plot(X{j}(k), j, 'ro', 'markersize', 2);
        set(gca, 'Xlim', [0 ddur]); % just show part of the train
    end;

    subplot(2,1,2);

    bar(hx, yh);
end;

return;

